Make prediction grid for relative prey weight models

Introduction

Make an evenly spaced UTM prediction grid with all spatially varying covariates for the stomach content data and biomass estimate data for the for diet and the biomass data

pkgs <- c("tidyverse", "tidylog", "sp", "raster", "devtools", "RCurl", "sdmTMB", "viridis", "terra", "ncdf4", "chron", "ncdf4", "tidync", "tidyterra") 

if(length(setdiff(pkgs, rownames(installed.packages()))) > 0){

    install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
  
}

invisible(lapply(pkgs, library, character.only = T))

# Source code for map plots
devtools::source_url("https://raw.githubusercontent.com/VThunell/Lammska_cod-fr/main/R/functions/map-plot.R")
# Packages not on CRAN
# devtools::install_github("seananderson/ggsidekick") # not on CRAN 
library(ggsidekick)
theme_set(theme_sleek())

# Set path
home <- here::here()

Read stomach data

d <- read_csv(paste0(home, "/data/stomach/stomachs_v1.csv"))

glimpse(d)

Rows: 132,274
Columns: 31
$ pred_ID            <chr> "LV_1968_9_19_1_37G4_52_NA_25", "LV_1973_5_11_15_37…
$ herring            <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ saduria            <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ sprat              <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ other              <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ other_inverts      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ other_chords       <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ rpw_tot            <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ rpw_tot_fb         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ rpw_sad            <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ rpw_spr            <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ rpw_her            <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ rpw_other          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ rpw_other_inverts  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ rpw_other_chords   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ year               <dbl> 1968, 1973, 1973, 1973, 1973, 1983, 1983, 1971, 197…
$ month              <dbl> 9, 5, 5, 5, 5, 12, 12, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6…
$ day                <dbl> 19, 11, 11, 11, 11, 16, 16, 19, 19, 19, 19, 19, 19,…
$ day_of_year        <dbl> 263, 131, 131, 131, 131, 350, 350, 170, 170, 170, 1…
$ pred_weight        <dbl> 129.63436, 38.25547, 38.25547, 45.19834, 32.05327, …
$ pred_weight_fb     <dbl> 152.97002, 46.30247, 46.30247, 55.27448, 38.37241, …
$ pred_length        <dbl> 25, 17, 17, 18, 16, 60, 77, 31, 33, 33, 30, 32, 28,…
$ pred_weight_source <chr> "estimated", "estimated", "estimated", "estimated",…
$ lat                <dbl> 54.25, 54.25, 54.25, 54.25, 54.25, 54.25, 54.25, 54…
$ lon                <dbl> 14.5, 15.5, 15.5, 15.5, 15.5, 15.5, 15.5, 14.5, 14.…
$ ices_rect          <chr> "37G4", "37G5", "37G5", "37G5", "37G5", "37G5", "37…
$ X                  <dbl> 467.4220, 532.5780, 532.5780, 532.5780, 532.5780, 5…
$ Y                  <dbl> 6011.453, 6011.453, 6011.453, 6011.453, 6011.453, 6…
$ data_source        <chr> "old_db", "old_db", "old_db", "old_db", "old_db", "…
$ Country            <chr> "LV", "LV", "LV", "LV", "LV", "LV", "LV", "LV", "LV…
$ Depth              <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…

Make the grid with depth

First make a grid for the biomass data, then subset that based on the extend of the stomach data

x <- d$X
y <- d$Y
z <- chull(x, y)

coords <- cbind(x[z], y[z])

coords <- rbind(coords, coords[1, ]) # close the loop

plot(coords[, 1] ~ coords[, 2]) # plot data

sp_poly <- sp::SpatialPolygons(
  list(sp::Polygons(list(sp::Polygon(coords)), ID = 1))
  )

sp_poly_df <- sp::SpatialPolygonsDataFrame(sp_poly,
                                           data = data.frame(ID = 1)
                                           )
cell_width <- 4 # increased from 3 to reduce prediction/index processing time and the NEMO hindcast has a spatial resolution of 3.7      

d |>
  summarise( n = n(), .by = month) # march is the most common month in the data. Use march for predictions

# A tibble: 12 × 2
   month     n
   <dbl> <int>
 1     9  7018
 2     5  4664
 3    12 11496
 4     6  9689
 5     4 25960
 6     3 35639
 7     1 10967
 8    11 15071
 9    10  2742
10     7  2193
11     2  6616
12     8   219

pred_grid_tmp <- expand.grid(
  X = seq(min(d$X), max(d$X), cell_width),
  Y = seq(min(d$Y), max(d$Y), cell_width),
  year = c(1963:2023)) |>
  mutate(month = 3)

ggplot(pred_grid_tmp |> filter(year == 2019), aes(X, Y)) +
  geom_point(size = 0.1) +
  theme_void() +
  coord_sf()

sp::coordinates(pred_grid_tmp) <- c("X", "Y")

inside <- !is.na(sp::over(pred_grid_tmp, as(sp_poly_df, "SpatialPolygons")))

pred_grid_tmp <- pred_grid_tmp[inside, ]

pred_grid_tmp <- as.data.frame(pred_grid_tmp)

ggplot(data = filter(pred_grid_tmp, year == 1999), aes(X*1000, Y*1000)) + 
  geom_point(size = 0.001, alpha = 0.5) +
  NULL

plot_map +
  geom_point(data = filter(pred_grid_tmp, year == 1999), aes(X*1000, Y*1000), size = 0.001, alpha = 0.5) +
  NULL

Warning: Removed 5091 rows containing missing values or values outside the scale range
(`geom_point()`).

# Add lat and lon
# Need to go from UTM to lat long for this one...
# https://stackoverflow.com/questions/30018098/how-to-convert-utm-coordinates-to-lat-and-long-in-r
xy <- as.matrix(pred_grid_tmp |> dplyr::select(X, Y) |> mutate(X = X*1000, Y = Y*1000))
v <- vect(xy, crs="+proj=utm +zone=33 +datum=WGS84  +units=m")
y <- project(v, "+proj=longlat +datum=WGS84")
lonlat <- geom(y)[, c("x", "y")]

pred_grid_tmp$lon <- lonlat[, 1]
pred_grid_tmp$lat <- lonlat[, 2]

ggplot(filter(pred_grid_tmp, year == 1999), aes(lon, lat)) + geom_point()

# Add depth now to remove islands and remaining land
# https://gis.stackexchange.com/questions/411261/read-multiple-layers-raster-from-ncdf-file-using-terra-package
# https://emodnet.ec.europa.eu/geoviewer/
dep_raster <- terra::rast(paste0(home, "/data/environment/Mean depth natural colour (with land).nc"))
class(dep_raster)

[1] "SpatRaster"
attr(,"package")
[1] "terra"

crs(dep_raster, proj = TRUE)

[1] "+proj=longlat +datum=WGS84 +no_defs"

plot(dep_raster)

pred_grid_tmp$depth <- terra::extract(dep_raster, pred_grid_tmp |> dplyr::select(lon, lat))$elevation

ggplot(pred_grid_tmp, aes(lon, lat, color = depth*-1)) + 
  geom_point()

pred_grid_tmp$depth <- pred_grid_tmp$depth*-1 # To make depth from negative elevation in relation to sea surface??

pred_grid_tmp <- pred_grid_tmp |> drop_na(depth)

pred_grid_tmp |> 
  filter(year == 1999) |> 
  drop_na(depth) |> 
  #mutate(water = ifelse(depth < 0.00000001, "N", "Y")) |> 
  ggplot(aes(X*1000, Y*1000, fill = depth)) + 
  geom_raster() +
  NULL

plot_map + 
  geom_point(data = pred_grid_tmp, aes(X*1000, Y*1000), size = 0.001) + 
  geom_sf() #Simple feature (https://r-spatial.github.io/sf/articles/sf1.html)

Warning: Removed 237961 rows containing missing values or values outside the scale range
(`geom_point()`).

plot_map + 
  geom_raster(data = filter(pred_grid_tmp, year == 1999), aes(X*1000, Y*1000, fill = depth)) + 
  geom_sf()

Warning: Removed 3910 rows containing missing values or values outside the scale range
(`geom_raster()`).

Combine oxy, temp, and sal from SMHI hindcast (and models of Copernicus and hindcast) for the pred_grid

# read Hindcast data
hindenv_df <- readRDS(file = paste0(home, "/data/environment/hindcast_1961_2017.rds")) |>
  filter(year > 1962) |>
  mutate(yearmonth = (year-1963)*12+month)

# We need to model the hindcast for year 2018-2023 and we use march as month.

# make a pred grid for those late observations
predhind <- hindenv_df |>
  distinct(lat, lon) |>
  replicate_df(time_name = "yearmonth", time_values = seq((2018-1963)*12+3,(2023-1963)*12+3, by=12)) |>
  mutate( model = as_factor("hindcast"),
          year = floor(1963+(yearmonth/12)),
          month = yearmonth %% 12) |> # mod, i.e. the remainder of an integer divide (here month)
  add_utm_columns(ll_names = c("lon", "lat"), utm_crs = 32633)

# Predict using the models for oxy, sal and temp. For reference see 1c_combine_oxysaltemp_data.qmd

# load model for oxygen
Mod_oxy <- readRDS(paste0(home, "/R/prepare-data/Mod_oxy.rds"))
# load model for salinity
Mod_sal <- readRDS(paste0(home, "/R/prepare-data/Mod_sal.rds"))
# load model for temperature
Mod_temp <- readRDS(paste0(home, "/R/prepare-data/Mod_temp.rds"))

oxypreds <- predict(Mod_oxy, newdata = predhind) |> mutate(oxy = est) |> dplyr::select(lat, lon, year, month, yearmonth, oxy)
salpreds <- predict(Mod_sal, newdata = predhind) |> mutate(sal = est) |> dplyr::select(lat, lon, year, month, yearmonth, sal)
temppreds <- predict(Mod_temp, newdata = predhind) |> mutate(temp = est) |> dplyr::select(lat, lon, year, month, yearmonth, temp)

# Combine hindcast preds for 2018-2023 with hindcast (1963-2017)
hindcast_allyears <- left_join(salpreds, oxypreds) |> left_join(temppreds) |> 
  bind_rows(hindenv_df |> dplyr::select(temp, sal, oxy, lat, lon, year, month, yearmonth))

Joining with `by = join_by(lat, lon, year, month, yearmonth)`
Joining with `by = join_by(lat, lon, year, month, yearmonth)`

map(hindcast_allyears, ~sum(is.na(.))) # No NAs

$lat
[1] 0

$lon
[1] 0

$year
[1] 0

$month
[1] 0

$yearmonth
[1] 0

$sal
[1] 0

$oxy
[1] 0

$temp
[1] 0

hindcast_allyears |> 
  summarise(nobs = n(), .by = c(month, year)) |>
  arrange(year)

# A tibble: 666 × 3
   month  year  nobs
   <dbl> <dbl> <int>
 1     1  1963 18300
 2     2  1963 18300
 3     3  1963 18300
 4     4  1963 18300
 5     5  1963 18300
 6     6  1963 18300
 7     7  1963 18300
 8     8  1963 18300
 9     9  1963 18300
10    10  1963 18300
# ℹ 656 more rows

hindcast_allyears |>
  summarise(mean_oxy = mean(oxy), .by = c(year, month)) |>
  ggplot() +
  geom_line(aes(year, mean_oxy, color = factor(month)), linetype = "dashed")

hindcast_allyears |>
  summarise(mean_temp = mean(temp), .by = c(year, month)) |>
  ggplot() +
  geom_line(aes(year, mean_temp, color = factor(month)), linetype = "dashed")

Extract oxy, sal and temp from hindcast_allyears to pred_grid

# add yearmonth
pred_grid_tmp2 <- pred_grid_tmp |>
  mutate(yearmonth = (year-1963)*12+month)

# function for extraction based on yearmonth
ext_envdat <- function(ayearmonth)  {

  ext_y = pred_grid_tmp2 |> filter(yearmonth == ayearmonth) |> dplyr::select(lon, lat) # coords from the stomach data, can only be lon and lat for extract()
  
  hindcast_allyears |>
    filter(yearmonth == ayearmonth) |>
    as_spatraster(xycols = c(2,1), crs = "WGS84", digits = 2) |>
    terra::extract(ext_y) |> 
    dplyr::select(oxy,sal,temp) |>
    bind_cols(pred_grid_tmp2 |> filter(yearmonth == ayearmonth)) 
    
}

# Use the ext_envdat function to extract env covs values to the observations for hindcast
pred_grid_tmp2 <- unique(pred_grid_tmp2 |> pull(yearmonth)) |> 
  map(\(x) ext_envdat(x)) |>
  list_rbind()

# Many observations lie just outside the hincast_allyears raster extent causing NA values of the covariates
theNAs <- pred_grid_tmp2 |> filter(is.na(oxy) | is.na(sal) | is.na(temp))
theNAs |> dplyr::select(lat,lon) |> distinct() |> summarise(n = n())# 427 latlon combos that are outside the hindcast extent

    n
1 427

# ... and they're along the coast which seems reasonable
plot_map_fc +
  geom_point(data = unique(theNAs |> dplyr::select(Y,X)), aes(X*1000, Y*1000)) +
  theme_sleek(base_size = 6) + 
  geom_sf() 

Warning: Removed 193 rows containing missing values or values outside the scale range
(`geom_point()`).

# We give them values based on the closest value in the hindcast 
# the 18300 unique hindcast latlon from which we can use get new lat lon to extract to predgrid with
hc_ull = hindcast_allyears |> filter(yearmonth == 1) |> dplyr::select(lon,lat)

# get the nearest neighbor
pl <- unique(theNAs |> dplyr::select(lon,lat) ) |>
  distance( y = hc_ull, lonlat = TRUE) 

# select which row in hc_ull has the smallest distance to each of theNAs
pl2 <- hc_ull[apply(pl, 1, FUN = which.min),] # cf a lengthy and slow dplyr version

# combine the real lat lons with their hc_ull corresponding neighbors, temporarily rename real lat lon for a left_join below.
pl3 <- bind_cols(unique(theNAs |> dplyr::select(lon,lat) |> rename(reallat = lat,
         reallon = lon)), pl2 )

# looks alright
plot_map_fc +
  geom_point(data = pl3 |> add_utm_columns(ll_names = c("reallon", "reallat"), utm_crs = 32633), aes(X*1000, Y*1000), color = "blue") +
  geom_point(data = pl3 |> add_utm_columns(ll_names = c("lon", "lat"), utm_crs = 32633), aes(X*1000, Y*1000), color = "red", alpha = 0.5) +
  theme_sleek(base_size = 6) + 
  geom_sf() 

Warning: Removed 193 rows containing missing values or values outside the scale range
(`geom_point()`).

Warning: Removed 194 rows containing missing values or values outside the scale range
(`geom_point()`).

# left join in the temporary pl3 latlons
pl4 <- left_join(theNAs |> rename(reallat = lat, reallon = lon), pl3)

Joining with `by = join_by(reallon, reallat)`

# left join in the covariates based on nearest neighbor lat lon and yearmonth
pl5 <- pl4 |> dplyr::select(-oxy,-sal,-temp,-year) |> 
  left_join(hindcast_allyears)

Joining with `by = join_by(month, yearmonth, lon, lat)`

pred_grid_tmp3 <- pl5 |> 
  dplyr::select(-lat,-lon) |> rename(lat = reallat, lon = reallon) |>
  bind_rows(pred_grid_tmp2 |> filter(!is.na(oxy) | !is.na(sal) | !is.na(temp)))

# all are there?(!)
pred_grid_tmp2 |> summarise(n()) == pred_grid_tmp2 |> summarise(n())

      n()
[1,] TRUE

# no NAs left
pred_grid_tmp3  |> filter(is.na(oxy) | is.na(sal) | is.na(temp))

 [1] X         Y         month     lon       lat       depth     yearmonth
 [8] year      sal       oxy       temp     
<0 rows> (or 0-length row.names)

# check plots
plot_map_fc +
  geom_tile(data = pred_grid_tmp3, aes(X*1000, Y*1000, fill = oxy), alpha = 0.5) +
  theme_sleek(base_size = 6) + 
  geom_sf() 

Warning: Removed 238510 rows containing missing values or values outside the scale range
(`geom_tile()`).

plot_map_fc +
  geom_tile(data = pred_grid_tmp3, aes(X*1000, Y*1000, fill = temp), alpha = 0.5) +
  theme_sleek(base_size = 6) + 
  geom_sf() 

Warning: Removed 238510 rows containing missing values or values outside the scale range
(`geom_tile()`).

plot_map_fc +
  #geom_tile(data = pred_grid_tmp3 |> filter(year %in% c(seq(1963,2023,10))), aes(X*1000, Y*1000, fill = sal)) +
  geom_tile(data = pred_grid_tmp3 |> filter(year == 2023), aes(X*1000, Y*1000, fill = sal)) +
  theme_sleek(base_size = 6) + 
  facet_wrap(~year) +
  scale_fill_viridis() +
  geom_sf() 

Warning: Removed 3910 rows containing missing values or values outside the scale range
(`geom_tile()`).

plot_map_fc +
  geom_tile(data = pred_grid_tmp3 |> filter(year %in% c(2015:2022)), aes(X*1000, Y*1000, fill = temp)) +
  theme_sleek(base_size = 6) + 
  facet_wrap(~year) +
  geom_sf() 

Warning: Removed 31280 rows containing missing values or values outside the scale range
(`geom_tile()`).

Add cod densities, only for post 1992!

# for the Mod, the fr data is instead used to scale the prediction grid
data_stats <- read_csv(paste0(home, "/data/survey/data_stats.csv")) 

Rows: 12328 Columns: 43
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (4): country, haul_id, ices_rect, month_year
dbl (39): year, lat, lon, quarter, month, sub_div, oxy, temp, sal, depth, X,...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

m3 <- readRDS(paste0(home, file = "/data/survey/m3.rds"))

# d_stats is used to scale the prediction grid when predictiong cod densities
data_stats_st <- data_stats |>
  summarise(depth_mean = mean(depth, na.omit = TRUE),
            depth_sd = sd(depth),
            sal_mean = mean(sal, na.omit = TRUE),
            sal_sd = sd(sal),
            oxy_mean = mean(oxy, na.omit = TRUE),
            oxy_sd = sd(oxy),
            temp_mean = mean(temp, na.omit = TRUE),
            temp_sd = sd(temp))

df_m3 <- pred_grid_tmp3 |>
  filter(between(year, 1993, 2021)) |>
  drop_na(oxy,
          sal,
          temp) |>
  mutate(depth = ifelse(depth < 0, 0, depth),
         depth_sc = (depth - data_stats_st$depth_mean)/data_stats_st$depth_sd,
         oxy_sc = (oxy - data_stats_st$oxy_mean)/data_stats_st$oxy_sd,
         sal_sc = (sal - data_stats_st$sal_mean)/data_stats_st$sal_sd,
         temp_sc = (temp - data_stats_st$temp_mean)/data_stats_st$temp_sd,
         depth_sq = depth_sc^2,
         temp_sq = temp_sc^2,
         year_f = as.factor(year),
         quarter_f = as.factor(1))

# Predict cod cpue density with model model from Max  
cpue_cod <- predict(m3, newdata = df_m3) |>
  mutate(log_density_cod = est) |>
  dplyr::select(lat, lon, month, year, log_density_cod) 

pred_grid_tmp4 <- pred_grid_tmp3 |>
  left_join(cpue_cod, by = c("month", "year", "lon", "lat"))

Add ICES areas

# https://stackoverflow.com/questions/34272309/extract-shapefile-value-to-point-with-r
# https://gis.ices.dk/sf/
shape <- shapefile(paste0(home, "/data/survey/ICES_StatRec_mapto_ICES_Areas/StatRec_map_Areas_Full_20170124.shp"))
head(shape)

  OBJECTID ID ICESNAME SOUTH WEST NORTH EAST AREA_KM2 stat_x stat_y Area_27
1        1 47     47A0  59.0  -44  59.5  -43     3178  -43.5  59.25  14.b.2
2        2 48     48A0  59.5  -44  60.0  -43     3132  -43.5  59.75  14.b.2
3        3 49     49A0  60.0  -44  60.5  -43     3085  -43.5  60.25  14.b.2
4        4 50     50A0  60.5  -44  61.0  -43     3038  -43.5  60.75  14.b.2
5        5 51     51A0  61.0  -44  61.5  -43     2991  -43.5  61.25  14.b.2
6        6 52     52A0  61.5  -44  62.0  -43     2944  -43.5  61.75  14.b.2
          Perc      MaxPer RNDMaxPer AreasList Shape_Leng Shape_Area
1 100.00000000 100.0000000       100    14.b.2          3        0.5
2  84.12674145  84.1267414        84    14.b.2          3        0.5
3  24.99803694  24.9980369        25    14.b.2          3        0.5
4  11.97744244  11.9774424        12    14.b.2          3        0.5
5   3.89717625   3.8971762         4    14.b.2          3        0.5
6   0.28578428   0.2857843         0    14.b.2          3        0.5

pts <- SpatialPoints(cbind(pred_grid_tmp3$lon, pred_grid_tmp3$lat), 
                     proj4string = CRS(proj4string(shape)))

pred_grid_tmp4$subdiv <- over(pts, shape)$Area_27

# Rename subdivisions to the more common names and do some more filtering (by sub div and area)
sort(unique(pred_grid_tmp4$subdiv))

[1] "3.d.24"   "3.d.25"   "3.d.26"   "3.d.27"   "3.d.28.1" "3.d.28.2" "3.d.29"  
[8] "3.d.32"  

pred_grid_tmp4 <- pred_grid_tmp4 |> 
  mutate(sub_div = factor(subdiv),
         sub_div = fct_recode(subdiv,
                              "24" = "3.d.24",
                              "25" = "3.d.25",
                              "26" = "3.d.26",
                              "27" = "3.d.27",
                              "28" = "3.d.28.1",
                              "28" = "3.d.28.2",
                              "29" = "3.d.29"),
         sub_div = as.character(sub_div)) |> 
  filter(sub_div %in% c("24", "25", "26", "27", "28", 2)) |> 
  filter(lat > 54 & lat < 59 & lon < 22)

# Add ICES rectangles
pred_grid_tmp4$ices_rect <- mapplots::ices.rect2(lon = pred_grid_tmp4$lon, lat = pred_grid_tmp4$lat)

plot_map +
  geom_raster(data = filter(pred_grid_tmp4, year == 1999), aes(X*1000, Y*1000, fill = depth)) +
  facet_wrap(~sub_div)

Warning: Removed 969 rows containing missing values or values outside the scale range
(`geom_raster()`).

pred_grid_done <- pred_grid_tmp4 |> dplyr::select(-subdiv)

Check

#zero valued env variables
zeropgplot <- pred_grid_done |>
  pivot_longer(cols=c(depth, oxy, sal, temp), names_to = "variable", values_to = "value" ) |>
  mutate(value_sc = scale(value), .by = variable) #|>
  #filter(value < 0)

plot_map +
  geom_raster(data = zeropgplot |> filter(variable == "oxy") |> filter(year %in% seq(1998,2023,5)), aes(X*1000, Y*1000, fill = value_sc)) +
  facet_wrap(~year) +
  scale_fill_viridis()

Warning: Removed 5814 rows containing missing values or values outside the scale range
(`geom_raster()`).

plot_map +
  geom_raster(data = zeropgplot |> filter(variable == "sal") |> filter(year %in% seq(1998,2023,5)), aes(X*1000, Y*1000, fill = value_sc)) +
  facet_wrap(~year) +
  scale_fill_viridis()

Warning: Removed 5814 rows containing missing values or values outside the scale range
(`geom_raster()`).

plot_map +
  geom_raster(data = zeropgplot |> filter(variable == "temp") |> filter(year %in% seq(1998,2023,5)), aes(X*1000, Y*1000, fill = value_sc)) +
  facet_wrap(~year) +
  scale_fill_viridis() # 1996 and 2006 was cold years

Warning: Removed 5814 rows containing missing values or values outside the scale range
(`geom_raster()`).

# cod density
plot_map +
  geom_tile(data = pred_grid_done |> filter(between(year, 1993,2021)) |> mutate(decade = round(year/10) * 10), aes(X*1000, Y*1000, fill = exp(log_density_cod))) +
  scale_fill_viridis(trans = "sqrt") +
facet_wrap(~decade)

Warning: Removed 28101 rows containing missing values or values outside the scale range
(`geom_tile()`).

Save

saveRDS(pred_grid_done, file = paste0(home, "/data/pred_grid.rds"))